Image registration is the process of determining correspondences between pixel elements in a pair of images that have common subject matter. It is an important technique in fields such as image matching, imaging device characterisation, and super-resolution. In image matching, two images are compared for common subject matter under the assumption that some geometrical transformation relates substantial portions of the two images.
Two images f1(x, y) and f2(x, y) can be related by a combination of Rotation, Scaling, and Translation (RST) transformations, such that:f2(x,y)=f1(s(x cos θ+y sin θ)+Δx,s(−x sin θ+y cos θ)+Δy)  (1)wherein s is a scale factor, θ is a rotation angle, and (Δx, Δy) are translations in x and y directions. The unknown rotation θ and scale s parameters may be determined from translation invariant representations Tf1 and Tf2 of the images f1(x, y) and f2(x, y), correspondingly. For example, a recombination of the Fourier magnitude |F| and the Laplacian of the Fourier phase ═2φ is invariant to any translation applied to the constituting image f:Tf=|F|+iA∇2φ  (2)where A is a scaling constant set to:A=max(|F|)/π  (3)to ensure that the recombined Fourier magnitude and phase information are roughly of equal magnitude.
After the rotation θ and scale s parameters are determined, the image f2(x, y) can be transformed to correct for the rotation and scaling. The translation parameters (Δx, Δy) can then be estimated from the image f1(x, y) and the transformed image f′2 (x, y) by finding a distinct peak in a cross-correlation image:C=ℑ−1(ℑ(f1)ℑ(f2)*)  (4)where ℑ(f) is the discrete Fourier transform of an image f(x, y), and ℑ(f2)* denotes the complex conjugation of the discrete Fourier transform ℑ(f2).
While a simple parametric transformation such as RST is suitable for registration of flat, rigid objects under frontal views, the transformation between two images in most applications is usually more complicated. In the general case, the transformation can be expressed as a free-form motion that warps one image onto a substantial part of the other image. This transformation is defined by a full warp map, which specifies a two-dimensional motion vector for each pixel of either image. This full warp map can be obtained from dense optical flow estimation. The full warp map can also be interpolated from a set of sparse point correspondences between the two images. Although there are many algorithms that detect and match salient points/regions from images, none of them is able to specify the level of registration accuracy that the algorithm returns.
Cramer-Rao Bound for Image Registration
The mean squared error of any estimate of a deterministic parameter in the presence of noise has a lower bound known as the Cramer-Rao Bound (CRB). Specifically, if a parameter vector m=[m0, m1, . . . mn]T is estimated from a given set of measurements, the CRB provides a lower bound on the error covariance matrix:
                              E          ⁡                      [                          ɛɛ              T                        ]                          ≥                                            (                              I                +                                                      ∂                    b                                                        ∂                    m                                                              )                        ⁢                                          F                                  -                  1                                            ⁡                              (                m                )                                      ⁢                                          (                                  I                  +                                                            ∂                      b                                                              ∂                      m                                                                      )                            T                                +                      bb            T                                              (        5        )            where ε={circumflex over (m)}−m is the estimation error (the hat sign denotes an estimate of the variable underneath); b=E[{circumflex over (m)}]−m is the bias of the estimator (E[.] is the expectation of the enclosed expression); I is the identity matrix; F (m) is the Fisher Information Matrix (FIM) that characterizes how well the unknown parameter vector m can be estimated from the observed data; and F−1 is the inverse of F. The ≧ sign in (5) means that the difference between the left matrix and the right matrix is non-negative definitive. As a result, the inequality holds for all diagonal terms.
If the estimator is unbiased (i.e., b=0), the expected variance of the parameters can be found directly from the main diagonal entries of the inverse matrix F−1:E[({circumflex over (m)}i−mi)2]≧[F−1(m)]ii  (6)
The Fisher information matrix is derived from the maximum likelihood principle. Let Pr(r|m) be the probability density function of an observed noisy data r(m), the Fisher information matrix is a measure of the steepness of the likelihood function around its peak:
                              F          ⁡                      [            m            ]                          ≥                                            E              [                                                ∂                                      ∂                    m                                                  ⁢                log                ⁢                                                                  ⁢                                  Pr                  ⁡                                      (                                          r                      |                      m                                        )                                                              ]                        [                                          ∂                                  ∂                  m                                            ⁢              log              ⁢                                                          ⁢                              Pr                ⁡                                  (                                      r                    |                    m                                    )                                                      ]                    T                                    (        7        )            
Since the peak of a steep likelihood function is less sensitive to noise than that of a smooth one, the FIM characterizes how precisely m can be estimated from the observed data.
Fisher Information Matrix for Image Registration
A direct image registration method searches for a parametric transformation between the coordinate systems of two images based on their intensity correlation. Assuming that both images I1 and I2 are noise corrupted versions of a noiseless scene I by two instances of zero-mean Gaussian noise with variance σn,I*1(x,y)=I1(x,y)+n1(x,y)=I(x,y)+n1(x,y)I*2(x,y)=I1(x,y)+n2(x,y)=I(x′,y′)+n2(x,y)  (8)where x′=f(x, y, m) and y′=g(x, y, m) are the coordinate transformations, and m=[m1, m2, . . . mn]T is the unknown registration parameter (e.g., under translation x′=x−tx, y′=y−ty and m=[tx, ty]T). Since the noise realizations n1 and n2 are normally distributed over the registration region S, the total probability of the unknown scene I given an estimate of m is:
                              Pr          ⁡                      (                          l              |                              m                ^                                      )                          =                              ∏            S                    ⁢                                          ⁢                                    1                                                σ                  n                                ⁢                                                      2                    ⁢                    π                                                                        ⁢                          (                              -                                                                            (                                                                        l                          1                          *                                                -                        l                                            )                                        2                                                        2                    ⁢                                                                                  ⁢                                          σ                      n                      2                                                                                  )                        ⁢                                          ∏                S                            ⁢                                                          ⁢                                                1                                                            σ                      n                                        ⁢                                                                  2                        ⁢                                                                                                  ⁢                        π                                                                                            ⁢                                  (                                      -                                                                                            (                                                                                    l                              2                              *                                                        -                                                          l                              ′                                                                                )                                                2                                                                    2                        ⁢                                                                                                  ⁢                                                  σ                          n                          2                                                                                                      )                                                                                        (        9        )            where the implicit coordinates for I1, I2 and I is (x, y) except for I′=I(x′, y′). The log-likelihood function therefore is:
                              log          ⁢                                          ⁢                      Pr            ⁡                          (                              l                |                                  m                  ^                                            )                                      =                                            -                              1                                  2                  ⁢                                                                          ⁢                                      σ                    n                    2                                                                        ⁢                                          ∑                S                            ⁢                              {                                                                            (                                                                        l                          1                          *                                                -                        l                                            )                                        2                                    +                                                            (                                                                        l                          2                          *                                                -                                                  l                          ′                                                                    )                                        2                                                  }                                              +          const                                    (        10        )            
From (7), the Fisher information matrix for a n-parameter vector m is thus a n×n matrix F with its entries computed as:
                                                                        F                ij                            =                            ⁢                                                -                                      E                    [                                                                                            ∂                          2                                                                                                      ∂                                                          m                              i                                                                                ⁢                                                      ∂                                                          m                              j                                                                                                                          ⁢                      log                      ⁢                                                                                          ⁢                                              Pr                        ⁡                                                  (                                                      l                            |                                                          m                              ^                                                                                )                                                                                      ]                                                  =                                  -                                      E                    [                                                                  ∂                                                  ∂                                                      m                            i                                                                                              ⁢                                              (                                                                              1                                                          σ                              n                              2                                                                                ⁢                                                                                    ∑                              S                                                        ⁢                                                                                          (                                                                                                      l                                    2                                                                    -                                                                      l                                    ′                                                                                                  )                                                            ⁢                                                                                                ∂                                                                      l                                    ′                                                                                                                                    ∂                                                                      m                                    j                                                                                                                                                                                                      )                                                              ]                                                                                                                          =                            ⁢                                                -                                      E                    [                                                                  1                                                  σ                          n                          2                                                                    ⁢                                                                        ∑                          S                                                ⁢                                                  (                                                                                                                    n                                2                                                            ⁢                                                                                                                                    ∂                                    2                                                                    ⁢                                                                      l                                    ′                                                                                                                                                                        ∂                                                                          m                                      i                                                                                                        ⁢                                                                      ∂                                                                          m                                      j                                                                                                                                                                                            -                                                                                                                            ∂                                                                      l                                    ′                                                                                                                                    ∂                                                                      m                                    i                                                                                                                              ⁢                                                                                                ∂                                                                      l                                    ′                                                                                                                                    ∂                                                                      m                                    j                                                                                                                                                                                )                                                                                      ]                                                  =                                                      1                                          σ                      n                      2                                                        ⁢                                                            ∑                      S                                        ⁢                                                                                            ∂                                                      l                            ′                                                                                                    ∂                                                      m                            i                                                                                              ⁢                                                                        ∂                                                      l                            ′                                                                                                    ∂                                                      m                            j                                                                                                                                                                                                      (        11        )            where the derivative of the noiseless image I′=I(x′, y′) with respect to each unknown parameter mi can be computed from its spatial derivatives and the registration model {x′, y′}={f(x, y, m), g(x, y, m)}:
                                          ∂                          l              ′                                            ∂                          m              i                                      =                                                            ∂                                  l                  ′                                                            ∂                                  x                  ′                                                      ⁢                                          ∂                                  x                  ′                                                            ∂                                  m                  i                                                              +                                                    ∂                                  l                  ′                                                            ∂                                  y                  ′                                                      ⁢                                          ∂                                  y                  ′                                                            ∂                                  m                  i                                                                                        (        12        )            Cramer-Rao Bound for 2D Shift Estimation
Using the general derivation of the Fisher information matrix in the previous subsection, the Cramer-Rao bound for any unbiased shift estimator can be derived. Two-dimensional (2D) shift estimation looks for a translational vector t=[tx, ty]T between the coordinate systems of the two images: x′=x−tx and y′=y−ty. The Fisher information matrix can be computed from (11) and (12):
                              F          ⁡                      (            t            )                          =                                            1                              σ                n                2                                      ⁡                          [                                                                                                                  ∑                        S                                            ⁢                                              l                        x                        2                                                                                                                                                ∑                        S                                            ⁢                                                                        l                          x                                                ⁢                                                  l                          y                                                                                                                                                                                                        ∑                        S                                            ⁢                                                                        l                          x                                                ⁢                                                  l                          y                                                                                                                                                                        ∑                        S                                            ⁢                                              l                        y                        2                                                                                                        ]                                =                                    1                              σ                n                2                                      ⁢            T                                              (        13        )            where lx=∂l′l∂x′=∂l′l∂x and ly=∂l′l∂y′=∂l′l∂y are spatial derivatives of the uncorrupted image I′. As can be seen in (13), the FIM for 2D shift estimation is proportional to a Gradient Structure Tensor T (GST) integrated over the region S.
Substitution of (13) into (6) yields the Cramer-Rao bounds of the variances of the registration parameters:
                                                        var              ⁡                              (                                  t                  x                                )                                      ≥                                          [                                  F                                      -                    1                                                  ]                            11                                =                                    σ              n              2                        ⁢                                          ∑                S                            ⁢                                                l                  y                  2                                /                                  Det                  ⁡                                      (                    T                    )                                                                                      ⁢                                  ⁢                                            var              ⁡                              (                                  t                  y                                )                                      ≥                                          [                                  F                                      -                    1                                                  ]                            22                                =                                    σ              n              2                        ⁢                                          ∑                S                            ⁢                                                l                  x                  2                                /                                  Det                  ⁡                                      (                    T                    )                                                                                                          (        14        )            where
      det    ⁡          (      T      )        =                    ∑        S            ⁢                        l          x          2                ⁢                              ∑            S                    ⁢                      l            y            2                                -                  (                              ∑                                          l                x                            ⁢                              l                y                                              S                )            2      is the determinant of T. Ignoring the second term of det(T), the Cramer-Rao bounds (14) are simplified to looser bounds:
                                          var            ⁡                          (                              t                x                            )                                ≥                                    σ              n              2                        /                                          ∑                S                            ⁢                              l                x                2                                                    ⁢                                  ⁢                              var            ⁡                          (                              t                y                            )                                ≥                                    σ              n              2                        /                                          ∑                S                            ⁢                              l                y                2                                                                        (        15        )            which clearly shows that the shift variance is linearly proportional to the input noise variance σn2 and inversely proportional to the total gradient energy in the shift direction. As a result, scenes with strong textures and little noise are likely to result in accurate shift estimation. However, the equality of the loose bound in (15) is hardly achievable (since
      (                  ∑        S            ⁢                        l          x                ⁢                  l          y                      )    2only vanishes when the orientation of maximum gradient energy is aligned with one of the grid axes).
Note that the CRB characterizes the shift variances based on an uncorrupted signal I, which is not available in practice. Fortunately, the total gradient energies of I can be approximated from those of the corrupted image I*1 and a noise instance n given a normal distribution N(0, σn):
                                          E            [                                          ∑                S                            ⁢                                                (                                      l                    x                    *                                    )                                2                                      ]                    =                                                                      ∑                  S                                ⁢                                  l                  x                  2                                            +                              E                [                                                      ∑                    S                                    ⁢                                      n                    x                    2                                                  ]                                      =                                                            ∑                  S                                ⁢                                  l                  x                  2                                            +                                                var                  ⁡                                      (                                          n                      x                                        )                                                  ⁢                                                      ∑                    S                                    ⁢                  1                                                                    ⁢                                  ⁢                              E            [                                          ∑                S                            ⁢                                                (                                      l                    y                    *                                    )                                2                                      ]                    =                                                                      ∑                  S                                ⁢                                  l                  y                  2                                            +                              E                [                                                      ∑                    S                                    ⁢                                      n                    y                    2                                                  ]                                      =                                                            ∑                  S                                ⁢                                  l                  y                  2                                            +                                                var                  ⁡                                      (                                          n                      y                                        )                                                  ⁢                                                      ∑                    S                                    ⁢                  1                                                                    ⁢                                  ⁢                              E            [                                          ∑                S                            ⁢                                                l                  x                  *                                ⁢                                  l                  y                  *                                                      ]                    =                                    ∑              S                        ⁢                                          l                x                            ⁢                              l                y                                                                        (        16        )            
where l*x=∂l*1 l∂x, l*y=∂l*1 l∂y, nx=∂nl∂x, and ny=∂nl∂y.
Cramer-Rao Bound for 2D Projective Registration
The Cramer-Rao bound is not only applicable to shift estimation, but also to more general motion models such as 2D affine and projective transformation. A 2D projective transformation, for example, is the motion of a static scene captured by a stationary camera or the motion of a moving planar surface. It is the most general planar motion model which includes translation, Euclidean, similarity, and affine transformations. Similarly to 2D translation, the Cramer-Rao bounds for the eight projective parameters are computed from an 8×8 Fisher information matrix.
Planar projective registration seeks an 8-parameter vector m=[m1, m2, . . . m8]T that transforms one coordinate system (x, y) into another (x′, y′):
                              [                                                    X                                                                    Y                                                                    D                                              ]                =                              [                                                                                m                    1                                                                                        m                    2                                                                                        m                    3                                                                                                                    m                    4                                                                                        m                    5                                                                                        m                    6                                                                                                                    m                    7                                                                                        m                    8                                                                    1                                                      ]                    ·                      [                                                            x                                                                              y                                                                              1                                                      ]                                              (        17        )                                                      x            ′                    =                                    X              D                        =                                                                                m                    1                                    ⁢                  x                                +                                                      m                    2                                    ⁢                  y                                +                                  m                  3                                                                                                  m                    7                                    ⁢                  x                                +                                                      m                    8                                    ⁢                  y                                +                1                                                    ,                              y            ′                    =                                    Y              D                        =                                                                                m                    4                                    ⁢                  x                                +                                                      m                    5                                    ⁢                  y                                +                                  m                  6                                                                                                  m                    7                                    ⁢                  x                                +                                                      m                    8                                    ⁢                  y                                +                1                                                                        (        18        )            Substituting
            [                                                                  ∂                                  x                  ′                                                            ∂                m                                                                                        ∂                                  y                  ′                                                            ∂                m                                                        ]        T    =            1      D        ⁡          [                                    x                                y                                1                                0                                0                                0                                                              -                                  x                  ′                                            ⁢              x                                                                          -                                  x                  ′                                            ⁢              y                                                            0                                0                                0                                x                                y                                1                                                              -                                  y                  ′                                            ⁢              x                                                                          -                                  y                  ′                                            ⁢              y                                          ]      into (12) yields:
                                          [                                          ∂                                  l                  ′                                                            ∂                m                                      ]                    T                =                              1            D                    [                                                                                          x                    ⁢                                                                  ∂                                                  l                          ′                                                                                            ∂                                                  x                          ′                                                                                                                                                          y                    ⁢                                                                  ∂                                                  l                          ′                                                                                            ∂                                                  y                          ′                                                                                                                                                                                ∂                                              l                        ′                                                                                    ∂                                              x                        ′                                                                                                                                  x                    ⁢                                                                  ∂                                                  l                          ′                                                                                            ∂                                                  y                          ′                                                                                                                                                          y                    ⁢                                                                  ∂                                                  l                          ′                                                                                            ∂                                                  y                          ′                                                                                                                                                                                ∂                                              l                        ′                                                                                    ∂                                              y                        ′                                                                                                              …                                                      ⁢                                                                                -                                          x                      ⁡                                              (                                                                                                            x                              ′                                                        ⁢                                                                                          ∂                                                                  l                                  ′                                                                                                                            ∂                                                                  x                                  ′                                                                                                                                              +                                                                                    y                              ′                                                        ⁢                                                                                          ∂                                                                  l                                  ′                                                                                                                            ∂                                                                  y                                  ′                                                                                                                                                                    )                                                                                                                                                        -                                              y                        ⁡                                                  (                                                                                                                    x                                ′                                                            ⁢                                                                                                ∂                                                                      l                                    ′                                                                                                                                    ∂                                                                      x                                    ′                                                                                                                                                        +                                                                                          y                                ′                                                            ⁢                                                                                                ∂                                                                      l                                    ′                                                                                                                                    ∂                                                                      y                                    ′                                                                                                                                                                                )                                                                                      ]                                                                                                          (        19        )            The 8×8 Fisher information matrix is rewritten from (11) and (19) as:
                              F          ⁡                      (            m            )                          =                              1                          σ              n              2                                ⁢                                    ∑              S                        ⁢                                                            [                                                            ∂                                              l                        ′                                                                                    ∂                      m                                                        ]                                ⁡                                  [                                                            ∂                                              l                        ′                                                                                    ∂                      m                                                        ]                                            T                                                          (        20        )            Due to a complex 8×8 matrix inversion, the exact formula for the Cramer-Rao bounds of 2D projective registration is not given here. The bounds can be computed from the diagonal entries of the inverse Fisher information matrix F−(m). Similarly to the shift estimation case, the lower variance bounds of the 2D projective parameters are proportional to the input noise variance and inversely proportional to the total gradient energy.
The Fisher information matrix, or alternatively the gradient structure tensor and hence the Cramer-Rao Lower Bound, quantifies the amount of information in an image for the determination of the n-parameter vector m. For shift estimation, the correlation information determines how precisely the shift can be estimated for a given area, given a certain amount of noise corrupting the image.